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Abstract. We introduce a novel technique for constructing higher-order variational integrators 
for Hamiltonian systems of ODEs. In particular, we are concerned with generating globally smooth 
approximations to solutions of a Hamiltonian system. Our construction of the discrete Lagrangian 
adopts Hermite interpolation polynomials and the Euler-Maclaurin quadrature formula, and in- 
volves applying collocation to the Euler-Lagrange equation and its prolongation. Considerable 
attention is devoted to the order analysis of the resulting variational integrators in terms of approx- 
imation properties of the Hermite polynomials and quadrature errors. A performance comparison 
is presented on a selection of these integrators. 



One of the major themes in geometric integration is symplectic methods for solving Hamiltonian 
systems of ordinary differential equations. Viewed as maps, such methods preserve the symplectic 
two-form underlying the dynamical evolution of the system. Variational integrators are an impor- 
tant class of symplectic integrators, which arise from discretizing Hamilton's principle. We refer 
the reader to [ll[ for a detailed discussion of the background theory of such methods. Variational 
integrators are automatically symplectic and momentum preserving. Moreover, they exhibit good 
energy behavior for exponentially long times. 

The construction of variational integrators combines the techniques from approximation theory 
and numerical quadrature and is related to the Galerkin approach of converting a differential 
operator equation into a discrete system. Galerkin methods are semi-analytic in the sense that the 
discrete solution is described by an element of a finite-dimensional function space, which provides 
an analytic expression for the numerical solution. Our goal is to develop variational integrators 
that would lead to globally smooth approximations of the solution. To pursue this goal, we adopt 
the space of piecewise Hermite polynomials in the Galerkin construction. It is a well-known result 
in approximation theory that Hermite polynomials allow for higher-order approximation of smooth 
functions at relatively low cost. Recall that Lagrange polynomials are a common choice in the 
construction of higher-order variational integrators. However, the computation of an interpolating 
polynomial in Lagrange form becomes unstable when the degree of the polynomial is high. Unlike 
Lagrange polynomials, higher-degree Hermite polynomials produce accurate and computationally 
stable results. 

Furthermore, the Galerkin construction based on piecewise Hermite polynomial interpolation 
allows us to adequately address the question of variational order error analysis. The order error 
analysis of variational integrators relies on determining the order with which a discrete Lagrangian 
L>d '■ Q x Q —> K approximates the exact discrete Lagrangian, 



where qoi(t) is a solution curve of the Euler-Lagrange equation that satisfies the boundary condi- 
tions goi(O) = <?0) Qoi(h) = qi- The standard way of performing the error analysis is by comparing 
the Taylor expansions of the exact discrete Lagrangian and the discrete Lagrangian. Consequently, 
the result critically depends on the extent to which the discrete trajectory is able to approximate 
the higher-derivatives of the exact Euler-Lagrange solution curve. To enable a robust variational 
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order error analysis for the Galerkin variational integrators, we propose a novel application of the 
collocation approach in the setting of discrete Lagrangian mechanics, which involves prolongations 
of the Euler-Lagrange vector field. The proposed approach has the advantage that one is able to 
prove optimal rates of convergence of the associated variational integrators as long as sufficiently 
accurate quadrature formulas are used. The proof crucially depends on the preliminary estimates 
on the higher-derivative approximation properties of prolongation-collocation curves. 

1.1. Outline of the Paper. We present a brief review of discrete variational mechanics and 
variational integrators in Section [21 and the Euler-Maclaurin quadrature formula in Section [3l In 
Section HI we introduce prolongation-collocation variational integrators, and perform variational 
order error analysis in Section [5j In Section [6l we present a few numerical examples, and present 
some conclusions and future directions in Section [7J 



2. Variational Integrators 

Let Q be the configuration manifold of a mechanical system, with generalized coordinates q. 
Consider the Lagrangian L : TQ — > R, where TQ is the tangent bundle of the configuration space 
Q. The tangent bundle has local coordinates (q,v). Further, let C(Q) = C([0, T],Q) denote the 
space of smooth trajectories q : [0, T] —> Q in the configuration manifold Q. The action integral 
S : C(Q) — >• R is defined as 

S(q) = ! T L(q(t),q(t))dt. 



JO 

The variational principle known as Hamilton's Principle states that 

5S = 0, for 8q(0) = 8q(T) = 0, 
which yields the Euler-Lagrange equations 

It is possible to rewrite the equations ([2]) in terms of the generalized coordinates and momenta (q,p) 
on the cotangent bundle T*Q (phase space). For this, we introduce the Legendre transformation 
FL : TQ -> T*Q, defined by 

The Hamiltonian H : T*Q — > R is given by 

H(q,p) = p-q-L(q,q)\ 3l . 

y dq 

One can show that equations ([2D are equivalent to Hamilton's equations (Theorem 1.3, p. 182 in 

B). 

dH . , dH , . 

(3) P = —Q-(p,q), q=-Q-(p,q). 

In the case of variational integrators, instead of discretizing the Euler-Lagrange equations ([2]), 
one discretizes Hamilton's principle. That is, one discretizes the action by introducing a discrete 
Lagrangian and replaces the action integral by an action sum, and applies the discrete version of 
the Hamilton's variational principle. The major advantage of this approach is that the resulting 
numerical algorithm automatically preserves the symplectic structure of the underlying dynamical 
system. 
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The discrete Lagrangian L d (qo,qx,h) is thought of as an approximation of the action integral 
along the curve segment between the points go ~ <?(0) and gi ~ q(h). Formally, this can be 
expressed as 

(4) L d (q , qi ,h) « / L(q(t),q(t))dt. 



We will neglect the /i-dependence and simply write L d (qo,qx) when it is not essential for the 
exposition of the material. Given the discrete sequence of times {t k = hk | k = 0, . . . , N}, h = T/N, 



a discrete curve in Q is denoted by {qk} k= Q, where g& ~ qifk)- The discrete action sum is a function 



that maps the discrete trajectories {qk}'k=o 



to 



N s 
k=0J 



I, and is given by 
N-l 

/, Ld(qk,qk+i)- 



k=0 



The discrete Hamilton's principle requires the discrete action to be stationary with respect to 
variations vanishing at k = and k = N. From this, one derives the discrete version of the 
Euler-Lagrange equations, which are known as the discrete Euler-Lagrange equations, 



(5) 



D 2 L d (q k - 1 , q k ) + DxL d (q k , q k+1 ) = 0, 



where k = 1,2, . . . ,N — 1. These equations implicitly define the one-step discrete Lagrangian map 
Fh d '■ Q x Q — > Q x Q- The discrete Legendre transforms F ± L ( ^ : Q x Q — > M are defined by 

F~L rf : (q ,qi) ^ (q , -D x L d {q Q , q x )), 
¥ + L d : (q ,qi) (qi,D 2 L d (q ,qi)), 

Pushing the discrete Lagrangian map Fi d forward to T*Q with the discrete Legendre transforms 
gives the discrete Hamiltonian map F^ d : T*Q — s> T*Q by Fj Jd = ¥ ± L d o o (¥ ± L d )^ 1 . 



fact that the definitions of Fl. are equivalent for the + and 
commutative diagram (Theorem 1.5.2 in [111]) 



The 

case is implied by the following 



(6) 



(<7o,£»o 




W2,P2j 



In coordinates, F^ d : (qo,Po) ^ (qi>Pi)> where 



(7) 



Po 



-DiL d (q , q{), px = D 2 L d {q ,q 1 ) 



A numerical quadrature can be used to approximate the integral in (JH). However, the functional 
form of the solution curve q(t) is required when applying a quadrature rule and it is, in general, 
unknown. In practice, one can choose an interpolating function on the interval [0, h] passing through 
go, qx- Then, a quadrature rule can be applied to the integral of the Lagrangian evaluated along the 
interpolating function. This approach fits the general framework of Galerkin integration methods. 
In more detail, for the construction of Galerkin Lagrangian variational integrators, one replaces the 
path space C([0,T],Q), which is an infinite-dimensional function space, with a finite-dimensional 
function space, C s ([0, T], Q). Commonly, one uses polynomial approximations to the trajectories, 
letting 

C s ([0, h],Q) = {q € C([0, h], Q) \ q is a polynomial of degree < s}. 
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IS 



An approximate action S(q) : C s ([0, h], Q) - 

s 

S{q) = h}^ biL(q(cih),q(cih)), 



i=l 



where Cj € [0,1] are quadrature points, 6j are quadrature weights, i = l,...,s. The Galerkin 
discrete Lagrangian is 



(8) 



Ld(qo,qi)= ext S(q). 

qeC s ([0,h],Q) 



In particular, for higher-order methods one takes q 6 C s ([0, h], Q) in the form 

s 

q(Th;q£,h) =^<?o^,s( r )> 



^=0 



where (?q, f = 1, ... ,s — 1, are the internal stages, I u ,s(t) are the Lagrange basis polynomials of 
degree s defined on the interval [0, 1]. Then, the integration scheme (qo,Po) >-> ((7i 5 Pi) is given by 



i=l 

s 

= /i fe, 

i=l 
s 

Pi = h 6j 



i=i 



^(cih)l 0>s (ci) + ^^{cih)l 0:S (a) 

9L, , , r , . 1 <9L ~ 
-Q-{cih)l^ a [ci) + -^-gT{Cih)l^ s (ci) 

^(cih)l s>s (ci) + ^^(cih)l SjS (a) 



, V = 1, s — 1 



It has been established (see [6|;lll|]) that the above scheme is equivalent to a symplectic partitioned 
Runge-Kutta method. Note that it yields a discrete solution that is, in general, only piecewise 
regular. In what follows, we introduce a construction of higher-order variational integrators with 
improved regularity across nodal times, which has as its natural variables the position, and its 
derivatives at only the nodal times, without the use of internal stages. 



3. Quadrature 

Here, we present a quadrature formula that we will use in our construction of a high-order discrete 
Lagrangian. The advantage of this particular rule is that it only involves function evaluations at 
the endpoints of the interval. When fast adaptive treecodes are used in conjunction with automatic 
differentiation techniques 15|], it is more efficient to obtain higher-order approximations using 
higher-derivative information at the endpoints, rather than evaluating the integrand at a number 
of internal stages. 

Theorem 1. (Euler Maclaurin quadrature formula) [1] If f is sufficiently differentiable 
(a, b), then for any m > 



Oil 



f(x)dx = - 



N-l 



f(a) + 2^2f(a + k9) + f(b) 



fc=l 



L>21 



,21 (/( 2 ^)(6)-/( 2 ^)(a)) 



Bo 



where B^ are the Bernoulli numbers, 6 = (6 — a)/N and £ S (a, b). 



>2m+2 jyg2m+3 i(2m+2) /v-\ 

(2m + 2)! J 
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Let us apply Theorem Q] to approximate an integral J Q f(x)dx in the simplest case when N = 1. 
It is easy to see that we obtain the following quadrature rule 



and the error of approximation is C(/i 2m+3 ). 

4. The Prolongation-Collocation Method 

In this section, we explain the construction of the discrete Lagrangian based on Hermite inter- 
polation and the Euler-Maclaurin quadrature formula. 

Motivation for Prolongation-Collocation Approach. The variational characterization of the 
exact discrete Lagrangian (pQ) naturally leads to the variational Galerkin discrete Lagrangian ([8]), 
where the infinite-dimensional function space C([0, h],Q) is replaced by a finite-dimensional sub- 
space, and the integral is approximated by a quadrature formula. While this leads to a com- 
putable discrete Lagrangian, one does not necessarily obtain an optimally accurate discrete La- 
grangian whose variational order is related to the best approximation properties of the chosen 
finite-dimensional function space. In particular, one finds that the variational Galerkin extremal 
curves do not necessarily approximate the higher-derivatives of the Euler-Lagrange solution curves 
with adequate accuracy. 

In retrospect, the fact that the variational Galerkin approach does not readily lead to computable 
discrete Lagrangians with provable approximation properties is not too surprising. By construction, 
variational Galerkin discrete Lagrangians associated with a sequence of finite-dimensional function 
spaces involve extremizers of a sequence of functionals. Since the sequence of finite-dimensional 
function spaces converges to C([0, h],Q), the sequence of functionals converges to the functional 
that appears in the variational characterization of the exact discrete Lagrangian. However, it is 
unclear that the sequence of extremizers converges to the extremizer of the limiting functional, since 
that corresponds to T-convergence [3] of the sequence of functionals. The issue of optimal rates of 
convergence of the computable discrete Lagrangians involves establishing rates of convergence of 
extremizers in terms of approximation rates of the finite-dimensional function spaces, which is an 
even more complicated process. 

As an alternative, we adopt the characterization of the exact discrete Lagrangian in terms of 
the Euler-Lagrange solution curve, and construct a discrete curve which approximates higher- 
derivatives of the Euler-Lagrange solution curve to an adequate level of accuracy. The latter is 
explored in detail in Section [5j 

Hermite Interpolation and Prolongation-Collocation. We commence by replacing q{t) in @ 
by its Hermite interpolant which is obtained by constructing a polynomial qd(t) such that values of 
q(t) and any number of its derivatives at given points are fitted by the corresponding function values 
and derivatives of qd(t)- In this paper we are concerned with fitting function values of q(t) and 
its derivatives at the end-points of the interval [0,h]. Consequently, a so-called two-point Hermite 
interpolant qd(t) of degree d = 2n — 1 can be used, which has the form 



(9) 



#(/) = £[/«>) + /(*)]-£ 



i=i 




n-l 



(10) 



Qd(t) 



Y, (q U) (0)H nj (t) + (-iyq^(h)H nj (h - tj) , 



3=0 



where 




n-j-l 




6 



MELVIN LEOK AND TATIANA SHINGEL 



are the Hermite basis functions. Note that for n = 1, the interpolant is a straight line joining q(0i) 
and q(h). By choosing one of the simple quadrature rules to discretize the integral in @ (e.g., the 
midpoint rule or trapezoidal rule), one obtains a class of well-known integrators which are at most 
second-order (see [III]). Therefore, the first nontrivial case of interest is n = 2, where we assume 
that the position and velocity data at the end points are available. From now on, we only consider 
n > 2 when applying the Hermite interpolation formula. The detailed derivation of (jlOp can be 
found, for example, in 0]. By construction, 

g W(0) = gW(0), q d r \h) = q^(h), r = 0, 1, . . . n - 1. 

Except for the step-size h, the discrete Lagrangian Ld(qo, qi, h) should only depend on qo « q(0), 
qi pa q(h). Therefore, letting q d (0) = qo and qd(h) = qi, we need to approximate the higher-order 
derivatives of q(t) by expressions that only depend on qo, q\. One natural approach, which is often 
found in the literature, is to use finite differences. In this work, we propose to apply the idea of 
collocation in conjunction with the Euler-Lagrange equations (|2|). The benefits of this approach 
will be exemplified later when discussing the variational error analysis of the proposed class of 
numerical integrators (see Section [5]) . 

The collocation approach [H is well-known in the theory of initial and boundary value problems 
for ODEs @;0j- Roughly speaking, the technique consists of determining the unknown parameters 
of a parameterized curve by requiring qd(t) to satisfy the ODE at a given set of points (collocation 
points). To define qd(t) uniquely, one sets the number of collocation points to be equal to the 
number of the available degrees of freedom. In our approach we use the method of collocation in a 
slightly unusual manner. In particular, since the parameters in (jlOp correspond to the derivatives 
of the solution curve q(t) at the end points of the interval [0, h], we are going to use t = and t = h 
as collocation points for the Euler-Lagrange equations (|2|) and consider the prolongation [13] of the 
Euler-Lagrange equations in order to generate a sufficient number of conditions. In other words, 
we increase the number of equations under consideration (not the number of collocation points) to 
match the number of degrees of freedom. 

For example, consider the case of the quintic Hermite interpolation, i.e., set n = 3 in (|10p . For 
separable Lagrangians of the form L(q,q) = \mq 2 — V(q), where m is the mass and V(q) is the 
potential energy term, the Euler-Lagrange equations © become a second-order ODE of the form 

m = /(?(<)), 

and its first-order prolongation can be expressed as 

q {3) (t) = f'(q(t))q(t). 

We set the boundary conditions q d (0) = qo an d q d (h) = Qi an d the collocation conditions 

&(0) = /(?d(0)), qf(0) = f'(q d (0))q d (0), 

q d (h) = f(q d (h)), q d 3 \h) = f \q d (h))q d (h). 

The above conditions constitute the system of six equations, which uniquely determines the fifth- 
degree polynomial q d (t) in the form of (JTUJ) . 

In general, for the Hermite polynomial of degree 2n— 1, one would need to differentiate the Euler- 
Lagrange equation n — 2 times, thus deriving a system of n — 1 equations for qd(t) , qy (t), . . . , q d n ^ (t) . 
Evaluated at and h, these (together with the Euler-Lagrange equations) give 2n — 2 collocation 
equations, which together with boundary conditions (q d (0) = qo,q d (h) = qi) constitute a sufficient 
number of conditions to determine the interpolant q d {t) uniquely. Note that for large n, the 
system of collocation conditions becomes nonlinear. Since the second and higher-order derivatives 
q d (0), q^\h) are given explicitly, the system can be recursively reduced to two implicit equations 
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involving <jd(0), Qd(h)- In this case, one would need to make use of a nonlinear root solver, such as 
the Newton-Raphson method, to determine qd(0),<jd(h). 

Further, in order to discretize the integral in @, we apply the Euler-Maclaurin quadrature 
formula Q. Recall that the formula involves derivatives of the integrand, in our case the Lagrangian 
L, with respect to the independent variable evaluated at the end-points of the integration interval. 
The latter, however, does not require the extensive computations typically associated with the 
interpolating polynomial due to the use of the Hermite interpolation formula and the collocation 
idea explained above. In more detail, we write 



L(q d (t),q d (t))dt « ~(L(q d (0),q d (Q)) + L(q d (h),q d (h))) 



m _ 

few* 



d 2l-l 

-jj^L(q d (t),q d (t)) 



t=h 



d 2l-l 

- d ^2TTL{Qd(t),qd(t)) 



t=0 



Provided that the degree of the interpolating polynomial is 2n — 1, we choose m < [n/2\, where 
the brackets denote the greatest integer lower bound for re/2. So for even n, \n/2\ = re/2 and for 
odd n, [n/2\ = (n — l)/2. This choice of m is justified by observing that the expressions for 

d 2l-l 



I = 1,2, ... , m, 



0,h, 



t=T 



include the derivatives of q d (t) up to order 2|_n/2j —1 + 1 < re, which satisfy the corresponding 
collocation conditions. 



Prolongation-Collocation Discrete Lagrangian. The Prolongation-Collocation discrete La- 
grangian is defined as follows, 

(12) L d (qo,qi,h) = ^(L(q d (0), q d (0)) + L(q d (h), q d (h))) 



\n/2\ 

ST fgj U21 

fe W 



j2l-l 



dt 21 



-L{q d (t),q d (t)) 



M-X 



t=h 



dt 21 



-L(q d (t),q d (t)) 



t=0 



where q d (t) € C S (Q) is determined by the boundary and prolongation-collocation conditions, 



(13) 



?d(0) = qo 
ffd(O) = /(go) 

i 3) (0) = /'(go)&(0) 



Qd(h) = qi, 
Qd(h) = /(gi), 

1 d 3 \h)=f(qi)q d (h), 



i \ d n 



t=o 



dt r - 



t=h 



One can include fewer than \n/2\ terms in the summation in (|12p . However, this will have an 
impact on the variational order of the corresponding integrator as will be further discussed in 
Section [5j The system of equations (fT3j) completely defines the discrete Lagrangian (fT2|) . Note 
that we used the second-order ODE q(t) = f(q(t)) as a prototype of the Euler-Lagrange equations 
for simplicity of notation only. The same idea applies for any smooth, not necessarily separable, 
Lagrangian function and the corresponding Euler-Lagrange equations. 

Given the initial conditions (qo,Po), the variational integrator has the form 

Pk = -DiL d (q k ,q k+1 ), 

(14) 

Pk+i = D 2 L d (q k ,qk+i), k = 0, 1, . . . , 
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and defines a one-step map (qk,Pk) ^ (Qk+i,Pk+i)- Generally, the equation pk = — -DiAife) <Zfc+i) 
together with the system of collocation conditions (|13[) can be reduced to a system of implicit equa- 
tions with respect to qk+i,Qd(tk), q<i(tk+i)- As soon as the solution is obtained via some appropriate 
nonlinear root-finding method, it is inserted into the equation pk+i = D2Ld(qk,qk+i)- 

When the Lagrangian has a relatively simple form, it makes sense to compute the expression for 
the discrete Lagrangian (|12p symbolically, which can be done using the symbolic module in Matlab 
or symbolic software such as Mathematica or Maple. Having computed a priori closed-form 
expressions for the right-hand side in (|14p , makes the implementation of the integrator particularly 
simple and fast. The reader is referred to Section [6] for some numerical examples. 

5. Variational Order Calculation 

The construction of variational integrators in the Galerkin framework naturally leads to the 
question of how it can be reconciled with the results from approximation theory of function spaces 
and numerical analysis of quadrature schemes. In particular, our goal is to explore the way the 
quantitative characteristics of the approximation errors enter the calculation of the convergence 
order of the respective integrators. Variational error analysis provides the right framework to 
pursue this goal. 



The variational error analysis introduced in and refined in |14l ]. is based on the idea that 



rather than considering how closely the numerical trajectory matches the exact flow, one can con- 
sider how the discrete Lagrangian approximates the exact discrete Lagrangian ([I]) which generates 
the exact flow map of the Euler-Lagrange equations. In other words, we are looking at the approx- 
imation error in 

L d (q(0),q(h))^ art / L(q(t),q(t))dt = Lf(q(0),q(h),h). 
q£C([0,h],Q) JO 

9(0)=90. 9C0=41 

We say that a given discrete Lagrangian is of order r if there exist an open subset U v C TQ with 
compact closure and constants C v and h v > so that 

(15) \\L d (q(0),q(h),h)-L$(q(0),q(h),h)\\ <C v h r+l 

for all solutions q(t) of the Euler-Lagrange equations with initial condition (q(0), q(Q)) £ U v and 
for all h <h v . In [11], the authors prove the equivalence of (|15p (cf. Theorem 2.3.1 in [111 ]) to: 

(i) the discrete Hamiltonian map Fi, d being of order r; 

(ii) the discrete Legendre transforms V^L d being of order r. 
In particular, the discrete Hamiltonian map is of order r if 

(16) \\F Ld (q(0)>P(P),h) -F Lf (q(0),p(Q),h)\\ < C v h r+ \ 

a 

for all solutions (q(t),p(t)) of the Hamilton's equations with initial condition (q(0),p(0)) € U w C 
T*Q and for all h < h w . The order of the discrete Legendre transforms is defined analogously. 
Recall from the diagram Q that 

Fi d ■ (<?o,Po) ^ {qi,Pi)- 

By construction, F l e (q(0) , q(h) , h) produces the values (q(h),p(h)) corresponding to the exact so- 

d 

lution of the Hamiltonian system ([3]), whereas Fi, d (q(0), q(h), h) produces the approximate values 
illi Pi) > ~ q{h), Pi ~ p(h). To summarize, the estimate (fTBI) provides the local order of con- 
vergence of the discrete trajectory (qk,Pk) to the exact flow (q(t),p(t)) of the Hamiltonian vector 
field. This order is the same as the order to which the discrete Lagrangian approximates the exact 
discrete Lagrangian, which we focus on. 

Before we explore the inequality (|15p for the Prolongation-Collocation discrete Lagrangian dis- 
cussed in Section HI we would like to establish the approximation error of q d \t) in comparison to 
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q^\t) for j = 1,2. . . , n, where q(t) is the exact solution of the Euler-Lagrange equation, and qd(t) 
is the Hermite interpolating polynomial (flC)]) of degree 2n — 1, constructed by letting q d (0) = q(0) 
and qd(h) = q{h), and imposing the prolongation-collocation conditions discussed in Section 0] at 
the endpoints. Note that this can be a difficult task in general, since the complexity of the col- 
location procedure escalates with the degree of the Hermite polynomial. However, we are only 
interested in the approximation order at the end-points of the interval [0,h], in which case the 
analysis is straightforward. 

Lemma 1. For q(t), qdif) as above, if qdij) = q(r) + 0(h p ) for some p > and r = 0,h, then 

q ( J ) (r) = q^(r)+0(hP), j = 3,...,n. 

Proof. Note that q(r) coincides with ^( T ) by construction. Therefore, we consider j starting 
from 3. As before, we restrict the proof to the case of a separable Lagrangian. Indeed, since the 
Euler-Lagrange equation is equivalent to the second-order ODE 

(17) q(t) = f(q(t)), 

it follows that for r = 0, h, 

q f {T) _ g (3) (r) = f'{ qd {r))q d (r) ~ f'(q(r))q(r) 

= f'(q(r))q d (r) - f(q(r))q(r) = f(q(r))(q d (r) - q(r)) = <D(hP\ 

provided there exists a uniform bound for /'. Consecutively differentiating (fTTl) and substituting 
the corresponding expressions for lower order derivatives, one can see that 9 W(t) (and qf {t)) can 
be represented as a polynomial in powers of q{r) (resp. qd{T~)) with coefficients which only depend 
on q(r) = q d {r), 

q (4) (r) = f"(q(T))q{rf + f'{q(r))f(q(r)) 

ff (5) (T) = /( 3) (g(r)M(r) 3 + g( T )( 3 /"(g(r))/(g(r)) + f'(q(r)) 2 ) 

As soon as / has bounded higher-order derivatives, the above formulas imply that the order of the 
approximation of q^\r) by q^\r) as a function of h is equal to the order of the approximation of 
q(r) by grf(r). We consider these expressions up to the j = n case, where n is determined by the 
number of collocation equations that were used to compute qd- □ 

In the next lemma, we determine the value of p in the relation 4rf(r) = q(r) + 0(h p ) in terms of 
the degree of the polynomial qd(t)- 

Lemma 2. Consider a polynomial qd(t) of degree d = 2n — 1 given by the formula 

n-l 

id(t) = K#nj(t) + (-iy aij H nd {h - 1)) , 

3=0 

where H n j (i) are the basis polynomial functions (|lip . By construction, 

qf (0) = a 0j , q ( j\h) = aij , j = 0, . . . n - 1. 

We let aoo = <z(0), aio = q(h). The coefficients a^j, a\j, j = 1,2, ...,n — 1 are obtained from 
the system of equations consisting of the Euler-Lagrange equation ([2]) and its prolongations. In 
particular, these are n — 1 differential equations evaluated on qd(t) at t = and t = h. Then for 
r = Q,h, 

qd(r) = q(r) + 0(h 2n - 1 ). 
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Proof. Let q(t) G C 2n ([0, h], Q) be the solution of the Euler-Lagrange equation ([2]) with boundary- 
conditions q(0) and q(h). Then, q(t) can be written in the form, 

q(t) = P 2n -3[q\(t) + Rn-Mtt), 



where 



n-2 



Pin-sm = E (? 0+2) (0)^-1,,^) + (-i)^ (j+2) (h)H n _ ld (h - *)) , 

i=o 

^-ifeK*) = . a)? ^ " ') • e G (0 ' 

and H n -i(t) are the basis polynomials defined by (jlip . See [3| for the proof of the above formula 
for sufficiently smooth functions. Note that q d (t) is a polynomial of degree 2n — 3 to which the 
same formula can be applied. In the latter case, the remainder term is identically zero. Hence, 

n-2 
3=0 

Subtracting q d (t) from q(t) gives 



*(*) = E (q^ +2 \0)Hn-i d (t) + (-iyq<j +2 \h)H n ^(h - t) 



(18) q(t)-q d (t) 

n-2 

= E ((^ ' +2) (°) - q { J +2 \0))H n -iAt) + (-iy(Q ij+2) (h) - q { J +2) (h))H n ^(h - t) 
3=1 

+ Rn-l[q\(t). 

Next, we integrate the above expression from to h. The left-hand side of the equation becomes 

f [q(t) - Ht)} dt = [q(h) - q d (h)} - [g(0) - q d (Q)] . 
Jo 

Further, observe that 

rh rh 

/ H n _ 1 j{t)dt= / H n - l>j (h-t)dt = C j h 2 , j = l,...n-2 
Jo Jo 

h 

t n-l( h _Qn-l = eh 271 " 1 , 



where Cj,C > are constants which do not depend on h. Now, let q d (0) = q(0) + 0(h p ) and 
q d (h) = q(h) + (D(h p ), where we wish to determine the value of p. It is easy to see that after 
integrating both sides of (|18p and rewriting it in terms of the order conditions we arrive at 

0{hF) = 0{h p+2 ) + o{h 2n - 1 ). 

Note that we used Lemma [T] to estimate the higher-order derivatives in (|18p . It follows immediately 
that p = 2n — 1 and the proof is finished. □ 

We are now ready to prove the main result of this section. 

Theorem 2. Assume that a Lagrangian function L : TQ — ^ M. is sufficiently smooth and its 
partial derivatives are uniformly bounded. Then, the discrete Lagrangian L d (qo,qi,h) constructed 
according to (|12p . (|13j) with qo = q(0) and q% = q(h), approximates the exact discrete Lagrangian 
Lj(q(0), q(h), h) with order 2\n/2\ + 3, for n > 3. In particular, when n is even, [n/2\ = n/2, 
and the order is n + 3. For odd n, \n/2\ = (n — 1) /2, and the order is n + 2. For n = 2, the order 
is equal to 3. 
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Proof. Observe that differentiability of the Lagrangian function L implies that it is Lipschitz con- 
tinuous in each of its arguments, given that the partial derivatives are uniformly bounded. We will 
make use of both, differentiability and Lipschitz continuity of L, in the proof below. 

We start with the simplest nontrivial case of n = 2, which corresponds to the space of piece- 
wise cubic polynomials in the Galerkin construction of a variational integrator. Note that it is 
sufficient to apply collocation to the Euler-Lagrange equation at the end-points of the interval 
[0,h] to uniquely define qd(t) satisfying the given boundary conditions. We obtain the order of 
approximation of the first derivative by applying Lemma [2 We use the simple trapezoidal rule 

h ^{f{a) + f(b)) + 0({b-af). 



f(t)dt 

to discretize the action integral in dU). Since L(q,q) is Lipschitz continuous, 
(19) L(q d (r),q d (r)) = L{q(r),q(j)) + 0(h 3 ). 

Therefore, 

h 



2 

£(£(g(0),5(0)) 



L d (q(0),q(h),h) = - (L(q d (0), q d (0)) + L(q d (h), q d (h))) 

L(q(h),q(h))) + 0(h A ) 

= L E d (q(Q),q(h))+0(h 3 ) + 0(h i ) 

= L*(q(Q),q(h))+0(h 3 ). 

The combination of the trapezoidal rule and the cubic Hermite interpolation makes the analysis of 
the approximation of the exact discrete Lagrangian by the discrete Lagrangian elementary. As we 
can see, the error in approximation is determined by the error of the quadrature rule. This is due 
to the fact that the derivatives are approximated to sufficiently high order. The pattern persists 
for higher-order Hermite interpolants due to our choice of quadrature method. 

It is straightforward to extend the above reasoning to the general case. Iteratively applying 
Lemma Q] and Lemma [2] to the derivatives of L d (q d (t),q d (t)) with respect to time t gives the 
relation 



d 21 



-i 



d 21 



-i 



t=0,h 



alt 21 



dt 2l^l L (ld(t),qd{t)) 
which is analogous to (fT9j) . Hence, 
L d (q(0),q(h),h) = ~(L(q d (0),q d (Q)) + L(q d (h),q d (h))) 



^L(q(t),q(t)) 



+ o(h 2n - 1 ), 



t=0,h 



K2J „ 

E — 

i=i 



21 



(21)1 



d 21 



-i 



dt 21 



^L(q d (t),q d (t)) 



d 21 



-1 



t=h 



dt 21 ' 



-L(q d (t),q d (t)) 



t=o 



^(L((0),q(0))+L( q (h), q (h)))+O(h 



2n\ 



K2J 

E 

l=i 



B21 h% 



(21)1 



d 2l ~ 



dt 21 



^L(q(t),q(t)) 



d 21 ' 



L%(q(0),q(h),h)+O(h 2 ^- 



t=h 
+ h.o.t. 



dt 2l ~ 



^L(q(t),q(t)) 



+ 0(h 2n ~ 1 ) 



t=0 



L^(q(0),q(h),h)+O(h 2 ^ 2 i +3 ). 



□ 



Remarks. Several remarks are in order. As we already noted, the error in approximation of 
the exact discrete Lagrangian by the Prolongation-Collocation discrete Lagrangian is determined 
by the order of the quadrature formula. In particular, the best order estimate is achieved if all 
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10 3 10~ 2 10 ' 10~ 4 10 3 10~ z 

time step cpu time 



Figure 1: Simple harmonic oscillator. Plots of: (a) the global errors for HEM and SRK4, (b) 
machine time versus the accuracy. The dotted line is the reference line for the exact order. 
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Figure 2: Simple harmonic oscillator. Energy error for the 4th order HEM and the two-stage 
symplectic RK4, step-size h = 0.2. 



the collocation conditions (I13j) are used. The collocation equations enter the terms under the 
summation in (|12p . which in turn determine to the order of accuracy of the quadrature formula. 

Secondly, if we write n = 2k, so that the degree of the interpolating polynomial is d = 2n — 1 = 
4k — 1, the choices n = 2k and n = 2k + 1 lead to the same maximal order of the quadrature, 
2k + 3. Therefore, it is preferable to use Hermite interpolating polynomials of order d = 4k — 1, 
which minimizes the computational effort for a discrete Lagrangian for a given order of accuracy. 

6. Examples 

6.1. Simple Harmonic Oscillator. We consider a harmonic oscillator system described by the 
equations 

q=p, p=-q. 

The total energy of the system is given by the Hamiltonian H(q,p) = ^p 2 + \q 2 ■ To test our 
method numerically, we used the 4th order variational integrator (HEM) constructed by means 
of the quintic Hermite interpolating polynomial and the Euler-Maclaurin quadrature formula. In 
the plots given in Figure [TJ Figure [H the resulting 4th order method is compared to the two-stage 
symplectic Runge-Kutta method of order 4. 
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Figure 3: Simple harmonic oscillator. The global approximation error of only the position trajec- 
tory. 




10 3 10 2 10 1 10 ' 10° 10' 

time step cpu time 



Figure 4: Planar pendulum. Plots of: (a) the global errors for HEM and SRK4, (b) machine time 
versus the accuracy. The dotted line is the reference line for the exact order. 



We would like to note the following interesting fact. Let us consider only the position component 
qk of the symplectic integrator (qk,Pk) ^ (<lk+i,Pk+i) an d compute the order of the corresponding 
one-step method qk >->■ qk+i- It turns out that in the case of the 4th order HEM method applied to 
the simple harmonic oscillator, the global error in the position component is of order 6. The error 
plots are given in Figure O where we can see that for the symplectic Runge-Kutta method the error 
remains to be of the 4th order. This does not however contradict the variational order analysis 
discussed in Section [5j The theorem mentioned therein establishes the order of the integrator 
ilkiPk) ^ (ik+iiPk+i) i n position-momentum variables, but allows the integrator (qk-i,qk) ^ 
{qkiQk+i) m position variables to have the same or higher order. 

6.2. Planar Pendulum. A planar pendulum of mass m = 1 with the massless rod of length / = 1 
is a Hamiltonian system for which the equations of motion are 

q = p, p = — sin q. 

In Figure HI we compare the performance of the method proposed in Section 0] with the two-stage 
symplectic Runge-Kutta method of the same order. HEM method encounters a slightly larger error 
in energy, but importantly it stays bounded for large time-intervals. 
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Figure 5: Pendulum. Energy error for the 4th order HEM and the two-stage symplectic RK4, 
step-size h = 0.2. 
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Figure 6: Duffing oscillator. Plots of: (a) the global errors for HEM and the Midpoint rule, (b) 
machine time versus the accuracy. The dotted line is the reference line for the exact order. 



6.3. Duffing Oscillator. The unforced undamped Duffing oscillator is a Hamiltonian system of 
equations 

q=p, p= q -q 3 , 

where the Hamiltonian function is 

H(q,p) = ^p 2 - ^q 2 + ^q 4 . 

The plots in Figure show the comparison of the 2nd order HEM variational integrator and the 
well-known Midpoint rule. The Midpoint rule is an implicit integrator and HEM is semi-implicit 
meaning that qt+i satisfies an implicit equation whereas Pk+i is computed explicitly. The plots 
demonstrate the superiority of the HEM over the Midpoint rule in terms of the computational time. 
The energy plots for both methods are given in Figure El 
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Figure 7: Duffing oscillator. Energy error for the 2nd order HEM and the Midpoint rule, step-size 
h = 0.2. 



7. Conclusions and Future Directions 

In this paper, we introduced a novel technique for constructing high-order variational integrators 
using collocation on the prolongation of the Euler-Lagrange equations. This relies on obtain- 
ing prolongation-collocation discrete curves with good approximation properties for the higher- 
derivatives. The resulting methods are particularly appropriate in combination with digital feed- 
back control, since they naturally yield position and its derivatives as the output, without the need 
to use interpolation in order to access such data. This also naturally leads to numerical trajectories 
with better regularity properties across the time nodes. 

It would be desirable to explore the connection between the methods proposed in this paper, 
and variational integrators based on global approximation techniques like splines, and to extend 
this work to the setting of Lie groups by incorporating techniques from Lie group variational 
integrators 0], and constructive approximation techniques on Lie groups Furthermore, it 



would be interesting to extend the prolongation-collocation techniques to Hamiltonian variational 
integrators [l(| and Hamilton-Pontryagin variational integrators 0] by considering prolongations 
of Hamilton's equations, and the implicit Euler-Lagrange equations. 
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